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VARIANCE  DETECTOR  REPORT 

INTRODUCTION 

This  report  documents  the  program   "vardet.c"  written   in   July,   1984   for 
post-   (FFT)   processor  activity/gap  detection.     The  algorithm  used  is  based  on 
the  assumption  that  the  variance  of  the  logs  of  the  frequency-  domain  magni- 
tudes of  a  signal    is  greater  than  the  corresponding  variance  of  the  channel 
noise. 

The  noise  is  assumed  to  be  quasi-ergodic,   that  is  ensemble  and  single  FFT 
averages  and  variances  are  assumed  to  be  the  same,  except  that  these  measures 
are  allowed  to  vary  slowly  in  time  using  a  weighted,   recursive  technique. 

Trade-offs  between  Type  I   (false  detection)  and  Type  II   (missed  signal) 
errors  can  be  controlled  with  an  interactively  defined  program  parameter. 

The  program  is  designed  to  warm-up  gracefully,   if  an  a  priori  estimate  of  an 
upper  bound  on  the  variance  of  the  logs  of  the  noise  magnitudes   is  available. 
In  the  absence  of  such  an  estimate,   the  warm-up  period  can  be  extended  inter- 
actively to  allow  the  program  to  learn  the  system  noise  levels. 

Additionally,   this  program  allows  the  addition  of  noise  to  an  input  signal 

so  that  program  parameters  may  be  tested  under  various  signal-to-noise  ratios. 

The  output  of  the  program  is   in   two  forms.      Indication  of  activity  or  gap  is 
printed  to  a   file  called   "GRAF"   for  graphing  on  the  line  printer.      For  periods 
in  which  activity  is  detected,   the  raw,   pre-FFT  data   is  printed   to  a   file 
called   "ACTIVE"   for  D/A  conversion   to  recover  the  gap-surpressed  signal.     An 
estimation  of  the  signal-to-noise  ratio   is  also  available  for  printing. 

TASK  DEFINITION 

The  problem  is  to  design  an  activity/ gap  detector  for  digitized  signal  data 
of  the  following  three  types:  speech,  fsk,  and  psk.  Running  time  is  limited 
by  the  specifications  of  the  mode  analysis  unit  into  which  the  detector  is  to 
be  incorporated,  and  number  of  operations  is  limited  by  hardware  considerations 
The  measure  of  effectiveness  can  be  given  by  the  complement  of  the  following 
function : 


FALSE  ALARMS  d(snr) 


1)  t  =   C  *       /  MISSED  DETECTIONS  d(snr)   +    / 

where  snr  is  the  signal-to-noise  ratio  and  the  limits  of  integration  are  over 
the  expected  operating  range,  probably  5  to  °°  dB.  This  equation  implies  that 
missed  detections  (Type  II  errors)  and  false  alarms  (Type  I  errors)  will  vary 
with  snr,  which  is  expected.  The  constant  C  implies  that  missed  detections 
and  false  alarms  might  be  weighted  differently  when  measuring  effectiveness. 
We  have  assumed  that  missed  detections  are  considerably  more  serious  than  false 
alarms. 

The  task  then,  in  summary,  is  to  develop  an  alogrithm  for  signal /gap  detection 
that  minimizes  0,  which  meeting  the  constraints  of  running  time  and  implementing 
hardware. 


THE  APPROACH 

It  was  noticed  upon  examination  of  typical  speech  and  fsk  data  that  the  variance 
of  the  logrithms  of  the  FFT  magnitudes  increases  when  signals  are  present  in  a 
channel.  A  detector  designed  to  exploit  this  finding  would  have  several  desir- 
able characteristics.  First,  variance  is  easy  to  compute,  requiring  only  the 
summing  then  squaring  and  squaring  then  summing  of  values  (about  3n  operations). 
Second,  the  variance  seems  to  be  a  fairly  sensitive  indicator  of  activity,  much 
more  sensitive  than  an  energy  level  (mean  value)  detector.  Third,  use  of  the 
logrithm  of  the  magnitude  of  the  FFT  makes  the  algorithm  independent  of  time- 
invariant  system  gain  levels.  Note  that: 

log(ax)  =  log  a  +  log  x  =  constant  +  log  x 
But  when  considering  the  variance, 

variance  [constant  +  f(x)]  =  variance  [f(x)] 

This  implies  that 

2)   variance  [log(ax)]  =  variance  [constant  +  log  x]  =  variance  [log  x] 

Thus  demonstrating  the  algorithm's  independence  of  any  constant  gain  a. 

The  use  of  logs  complicates  the  algorithm,  increasing  its  running  time  and  pre- 
venting the  calculation  of  the  measure  variance  [logl I FFT(signal  +  noise)  I  I D 
as  a  function  of  the  signal  and  noise  distribution  parameters. 

THE  ALGORITHM 

The  variance-based  detection  algorithm  can   be  given  as: 

If  variance  [log    I  I FFT(waveform) | | ]   >  some  threshold 
declare  activity. 

Otherwise,   declare   inactivity. 

The  problem  becomes  establishing  an  appropriate  threshold.      Clearly,   a   low 

threshold  will   lead  to  a   high  false  alarm  rate,  while  a   high  threshold  will 

cause  missed  detections.     The  threshold  must  be  related   to  the  variance  of 
the  noise  alone  as   follows: 

threshold   >  variance  [log    I  I FFT(noise) I  I ] 

Hereafter,   the  right  hand  side  above  will    be  referred  to  as  var(noise)   for 
simpl  icity. 

The  communication  channels  under  consideration  are  4   kHz  wide,   but  only  the 
300  -  3700  Hz  band  contains  activity.     This  means  that  the   first  and  last 
300  Hz  bands   in   the  channel    should  contain  only  noise.      So  with  a   64  point 
FFT,   the  middle  56  bins  will   contain  the  activity,   if  present  and  the   first 
and  last  4  bins  will    contain   the  FFT  of  the  noise  only.      In   practice,  we  have 


noticed  that  the  two  bins  immediately  adjacent  to  the  activity  band  on  each 
side  are  contaminated  due  to  filter  and  windowing  imprecision.  This  leaves  a 
total  of  4  bins,  2  above  and  2  below  the  activity  channel,  as  measures  of  the 
noise  only. 

Four  samples,  however,  are  too  few  to  get  an  adequate  measure  of  the  noise 
variance.  If  we  can  assume  the  noise  to  be  ergodic,  we  can  use  the  variance 
over  time  as  a  measure  of  the  variance  over  frequency.  That  is,  if  the  noise 
is  ergodic,  the  statistical  measures  at  any  time  over  all  frequencies  will  be 
the  same  as  the  measures  at  any  frequency  over  all  time.  Because  we  have  4 
noise  bins,  we  will  have  4  independent  measures  of  noise  variance  over  time. 
Given  the  assumption  of  ergodicity,  these  will  be  close  to  the  variance  of  the 
noise  over  frequency  and  can  be  used  in  setting  the  threshold. 

We  experimented  with  different  ways  of  using  these  4  measures  to  arrive  at  a 
single  estimate  of  var(noise).  One  method  used,  but  discarded,  is  to  take  the 
maximum  of  the  4  values  at  any  given  time.  The  problem  with  this  method  is 
the  extreme  variability  between  successive  threshold  values.  A  second  method, 
found  to  be  more  stable,  is  to  take  the  average  of  these  4  variances  as  an 
estimate  of  the  true  variance.  To  establish  a  threshold,  we  increase  this 
value  by  multiplying  it  by  a  "kluge"  factor  (K)  to  be  determined  empirically 
and  depending  upon  the  distribution  of  the  noise  over  time  and  the  trade-off 
between  Type  I  and  Type  II  errors.  Thus, 

3)   threshold  =  K  *  average  over  4  bins  [var(noise)  over  time]  >  var(noise) 

We  can  summarize  the  detection  algorithm  as: 

If  variance  across  inner  56  activity  bins  at  any  time  [log  I  I FFT(waveform) I  I ] 

is  greater  than 

K  *  average  across  the  4  outer  bins  {variance  over  time  of  each  of  the  4  bins 
[log  I  I  FFT(waveform)  I  |  ]  } 

declare  activity, 

Note  again  that  the  variance  across  the  inner  activity  bins  is  at  a  single  time 
window,  while  the  variance  of  the  noise  is  measured  across  all  time  windows. 

In  application,  we  are  using  a  64  point  FFT  requiring  128  data  points.  In 
general,  activity  will  not  truely  begin  or  end  at  the  start  of  one  of  these  FFT 
cycles.  Therefore,  we  expect  a  one  cycle  uncertainty  regarding  the  true  start 
and  end  of  activity.  A  data  sampling  rate  of  8000  for  the  4000  Hz  transmission 
channel  requires  62.5  FFT  cycles  per  second  of  128  data  points  for  real  time 
processing.  Thus,  one  cycle  would  translate  to  a  1 6  millisecond  uncertainty 
regarding  the  true  start  and  end  of  activity.  Speech  activity  often  begins  and 
ends  at  low  energy  and,  presumably,  low  variance  levels,  so  the  true  uncertainty 
in  practice  would  be  even  higher. 

ALLOWING  FOR  TIME-VARYING  NOISE  DISTRIBUTIONS 

We  would  like  to  make  as  few  assumptions  regarding  the  noise  distribution  as 
possible.  We  have  thus  far  assumed  only  ergodicity.   It  may  be  possible  to 


create  an  algorithm  which  is  robust  against  violations  of  this  assumption  by 
allowing  for  a  slowly  (with  respect  to  the  time  of  one  FFT)  time-varying  noise 
distribution.  This  could  be  accomplished  by  using  a  moving,  weighted  measure 
of  the  variance  across  time  of  each  of  the  4  noise  bins. 

We  will  start  by  finding  the  moving,  weighted  mean  of  the  magnitude  of  the  FFT 
for  each  of  the  4  noise  bins. 

n     "a°  +~  a"""a"". "...."..  an 

where  x  is  the  magnitude  of  the  FFT  of  the  i  th  bin  at  each  time  window  and 
o  <  a  <  1  is  the  attenuation  constant. 

For  the  purposes  of  calculation,  this  can  be  written  recursively  as 

a. 
^        a*X    +   x 
5)       X  n n_ 

n+1      d  n+1 


d     =   a*d  +  1 
n+1         n 


Now  we  calculate  the  moving,  weighted  variance  based  on  this  weighted  mean, 


6)   VAR  = 


a°  (xn  -  Xn)2  +  ai(x    _  xn_2)2  + a"  (x0  -  X0) 


2 


a0  +  a1  +  a2  +  a 


n 


This  can  be  written  recursively  as 

a, 

7)       U  =     a*U     +     (x  -   X       )2 

n+1  n  n+1  n+1 

d  =     a*d       +     1 

n+1  n 

U 

VAR         =     n+1 


n+1  (n+l)*d 

n+1 

The  use  of  a  moving,  weighted  variance,   however,   has  effects  not  anticipated 


Basing  the  variance  on  a  moving  mean  biases  the  measure,  causing  it  to  be  an 
underestimate.  There  is  a  mechanical  analogy  illustrating  this  effect.  If  a 
pendulum  is  allowed  to  swing  on  a  moving  pivot,  the  pendulum  will  swing  faster. 
The  pivot  moves  in  the  direction  of  the  center  of  gravity  (moving  mean).  The 
faster  swing  indicates  that  the  momemt  of  angular  momentum  (second  moment  or 
variance)  has  been  lowered. 

Even  if  the  noise  distribution  were  completely  specified,  finding  the  degree 
to  which  the  variance  is  underestimated  by  this  method  would  still  be  a  for- 
midable task.  Recall  from  equation  3),  however,  that  the  average  of  the  4  bin 
variances  will  be  multiplied  by  the  kluge  factor  K.  If  the  noise  distribution 
doesn't  vary  too  wildly,  the  kluge  factor  can  be  turned  empirically  to  com- 
pensate for  the  underestimate  of  the  variance. 

It  will  be  noted  that  the  denominator  in  both  the  moving  mean  and  variance  cal- 
culations is  a  geometric  series  with  a  constant  0  <  a  <  1.  This  series  converges 
to  l/(l-a).  For  the  purposes  of  implementing  this  algorithm  into  hardware,  this 
limiting  value  is  used  through  all  calculations.  This,  of  course,  means  that 
for  the  early  cycles  the  mean  will  be  grossly  underestimated  and  the  variance 
will  be  grossly  overestimated. 

GRACEFUL  WARMUP 

It  was  mentioned  above  that  approximating  the  denominator  of  the  moving  statis- 
tics by  the  limit  of  the  geometric  series  causes  early  measures  to  be  yery   poor. 
There  is  a  second  reason  to  expect  problems  in  the  early  cycles. 

a. 

Note  that  in  equation  5)  the  last  term  on  the  right  hand  side  (x0  -  X0)  is 
always  zero  because 

x0 
X  =  - 

0    1 

This  means  that  VAR  =  0  and,  by  equation  3),  the  threshold  =  0.   For  small 
sample  sizes  (shortly  after  the  algorithm  is  turned  on)  the  variance,  and  hence 
the  threshold,  will  be  grossly  underestimated.  Consequently,  early  waveforms 
will  invariably  declared  activity,  whether  or  not  activity  is  truely  present. 

This  is  compensated  for  in  two  ways.   First,  the  algorithm  can  be  given  an 
adequate  chance  to  warm  up.   In  the  computer  program  to  be  presented  in  following 
sections,  the  warm  up  period  is  interactively  specified  by  the  user.  Second, 
if  an  upper  bound  on  var(noise)  is  available  a  priori,  it  can  be  used  to  replace 
(x  -  X)  in  equation  6).  The  following  program  also  incorporates  this  method. 

NOISE  TRANSIENT  AND  ACTIVITY  DROP-OUT  FILTERING 

We  noticed  in  application  that  in  the  midst  of  a  string  of  non-activity  de- 
clarations, a  spurious  declaration  of  activity  would  often  take  place.   In 
actuality,  speech  activity  of  single  cycle  duration  would  be  most  unusual. 
Similarly,  single  cycle  drop-outs  would  often  occur  in  the  midst  of  a  string 
of  activity.  To  eliminate  these  problems,  we  adopted  a  "2  out  of  3"  filtering 
criterion.  The  filter  is  activated  at  the  discretion  of  the  program  user  by 


interactively  responding  to  query  regarding  filtering.     The  filter  simply  allows 
activity  to  be  declared  only  if  2  of  the  last  3  cycles  indicated  activity. 
Similarly,  non-activity  is  declared  only  if  2  of  the  last  3  cycles   indicated 
non-activity. 

One  effect  of  the  filter  is  to  cause  the  first  cycle  of  a  string  of  activity  to 
be  invariably  declared  as  non-activity.      If  real-time  delay  and  storage  require- 
ments are  not  a  problem,  all   data  may  be  saved  through  the  following  cycle.      If 
a  cycle  is  the  first  in  a  string,   the  data  may  be  retrieved  after  the  next  cycle 
and  then  written  to  the  output  files. 

ESTABLISHING  OPTIMAL  PARAMETER  VALUES 

Parameters  K,  a,  and  the  initial  value  of  the  var(noise)  can  be  optimally  set 
by  exhaustion  using  equation  1).  Due  to  the  magnitude  of  this  project,  it  has 
not  yet  been  done.  Rather,  preliminary  values  of  these  parameters  have  been 
found  which  yield  satisfactory  results  for  the  waveforms  thus  far  tested. 

Under  the  previous  assumptions,  the  value  of  K  should  be  greater  than  1. 
Experimental  evidence,  however,  has  lead  us  to  believe  that  for  signal  and  noise 
environments  thus  far  tested,  K  should  be  less  than  1  for  optimal  performance. 
There  are  two  reasons  for  this  surprising  result.  The  first  is  the  judgment 
that  Type  II  errors  (missed  detections)  are  more  serious  than  Type  I  errors 
(false  alarms).  Recall  that  lowering  the  threshold  decreases  missed  detections 
at  the  cost  of  false  alarms.  The  second  reason,  however,  is  fundamentally 
related  to  the  characteristics  of  the  noise.  We  noticed  that  when  activity  is 
present,  the  magnitude  of  FFT's  of  the  4  outer  bins  drops  about  40%  for  the 
waveforms  tested.  This  would  be  true  if  an  automatic  gain  control  were  present 
somewhere  in  the  system.  A  time-varying  gain  causes  the  noise  to  be  non-ergodic. 
The  variances  over  time  of  the  outer  bin  FFT's  increase  because  the  magnitudes 
change  (decrease)  whenever  activity  is  present.  All  this  serves  to  overestimate 
the  noise  variance.  Thus,  it  appears  that  an  optimal  value  of  K  will  be  some- 
what below  1.   In  the  following  program,  the  user  specifies  K  interactively. 

The  attenuation  constant  must  also  be  specified,  by  nature,  between  0  and  1. 
A  low  value  allows  greater  movement  of  the  weighted  mean,  thus  lowering  the 
measure  of  the  variance  and  causing  a  high  false  alarm  rate,  which  can  be  com- 
pensated for  by  raising  the  value  of  K.  A  high  value  causes  the  mean  and  variance 
to  become  more  insensitive  to  changes  in  the  noise  parameters.  This  may  be 
desirable  as  it  causes  the  threshold  to  remain  stable  through  brief  periods  of 
noise  instability. 

The  parameters  were  chosen  on  the  basis  of  tests  on  a  102400  value  data  set. 
This  set  yielded  800  FFT  cycles.  Approximately  258  of  these  cycles  were 
identified  ahead  of  time  as  containing  activity.  The  remaining  552  were  felt 
to  be  non-active.  We  found  that  for  the  single  data  set  tested,  with  no  added 
noise,  the  values  K  =  0.8,  a  =  0.997  and  U  =  57.0  worked  well,  yielding  roughly 
41  missed  detections  and  26  false  alarms.  For  the  purposes  of  speech  compres- 
sion, this  result  is  entirely  satisfactory. 

THE  DETECTION  PROGRAM 

A  flow  chart  of  the  program  is  given  as  Figure  1.  The  program  itself,  written 
in  the  C  programming  language,  is  included  in  the  appendix  along  with  a  variable 


list.  The  main  program  requires  four  external  functions  (subroutines).  These 
are:  iread  and  iwrite,  both  written  by  Steven  D.  Beck  to  read  into  the  main 
program  the  digitized  signal  values  and  to  write  to  an  output  file  the  gap  sur- 
passed digitized  signal;  fft2,  provided  by  Steven  Beck,  to  perform  the  Fast 
Fourier  Transform;  and  actgraf2,  written  by  Walter  Underwood  to  create  a  file 
for  graphing  on  the  line  printer  the  indication  of  gap  or  activity. 

NOISE  SIMULATION 

The  program  has  the  capability  of  adding  noise  to  the  signal  prior  to  the 
variance  calculations  so  that  program  performance  may  be  tested  at  low  signal- 
to-noise  ratios.  The  program  assumes  that  the  warmup  period  consists  of  non- 
activity.  The  complex  FFT  data  for  each  warmup  cycle  is  stored.  We  exploit  the 
linearity  of  Fourier  transforms  to  add  the  warmup  FFTs  multiplied  by  an  inter- 
actively set  gain  to  the  FFTs  of  subsequent  cycles  to  get  the  equivalent  FFTs 
of  noisy  data.  The  warmup  FFTs  are  added  to  the  subsequent  FFTs  cyclically. 

The  advantage  of  this  method  is  that  no  assumptions  regarding  noise,  except 
stationarity,  are  made.  We  are  assuming,  however,  that  the  warmup  period 
contains  noise  only.  Consequently,  the  noise  simulation  works  only  with  data 
files  that  begin  with  non-activity. 

ALGORITHM  PERFORMANCE 

As  previously  mentioned,  parameter  values  were  found  which  yielded  approximately 
a  16%  missed  detection  rate  and  a  10%  false  alarm  rate,  basing  both  percentages 
on  the  number  of  cycles  actually  felt  to  contain  activity.  These  values  were 
found  with  data  exhibiting  approximately  a  20  dB  signal-to-noise  ratio  for  single 
activity  cycles.  Added  noise  caused  graceful  performance  degradation,  primarily 
increasing  the  false  alarm  rate. 

An  obvious  weakness  of  the  algorithm,  however,  was  its  inability  to  detect  the 
spoken  word  "six"  as  activity.  The  "s"  and  "x"  sibilate  sounds  are  broadband. 
The  effect  of  these  sounds  is  identical  to  adding  noise  to  the  activity  band. 
The  variance  of  the  FFT  of  this  band  is  not  increased  by  the  presence  of  these 
sounds,  hence  their  undetectabil  ity.  This  appears  to  be  a  fundamental  weakness 
of  the  algorithm  and  not  the  result  of  poorly  chosen  detection  parameters. 

HARDWARE  IMPLEMENTATION 

A  flow  chart  indicating  possible  hardware  implementation  is  attached. 

SUMMARY 

A  detection  algorithm  has   been  written  that  gives   good  results   for  speech 
activity  and  also   is   simple  enough  for  hardware  implementation.      Performance 
with   fsk  and   psk  signals   has  not  yet  been  tested,   but  we  expect  satisfactory 
performance  with  these  signal    types  as  well.      No  assumptions   regarding  the 
distribution  of  the  underlying  noise  are  made,   although  alogrithm  parameters 
are  tuned  on  the  basis  of  empirical    knowledge  of  the  noise.     The  system  is   ro- 
bust against  changing  noise  parameters  and  degrades   gracefully  with  decreasing 
signal-to-noise  ratio.     The  algorithm  has  proven  unable  to  detect  sibilates, 
but  otherwise  gives  excellent  performance. 


APPENDIX 


The  following  appendix  contains  the  flow  chart,  the  variable  list  and  a  listing 
of  the  program  "vardet.c". 


VARIABLE  LIST 


ang 
a  ten 

avesnr 

avevar 
counter 


countO 
countl 
count2 


denom 

df 

end 
fin 

fp 

fph 

fpw 

h[128] 

i 

idf 
im 

inoise[j][64] 

irec 
isactive 


added  noise  gain,  used  when  simulating  added  noise 

atenuation  constant  for  moving  weighted  mean  and 
variance.  0  <  aten  <  1 

average  signal-to-noise  ratio  over  6  declarations  of 
activity 

average  of  var(noise)  across  the  s4  outer  bins 

index  used  when  windowing  digitized  waveform  prior 
to  FFT 

flags  used   for  2  out  of  3   filtering.      l=activity, 
0=inactivity.      CountO  is  the   flag  for  the  most  recent 
FFT,   countl    for  the  previous,   count2   for  second 
previous. 

denominator  in  weighted,  moving  mean  and  variance. 
Geometric  series  approximated  by  limit  value  1/(1 -a). 

degrees  of  freedom  of  inner  bin  variance  =  #  of 
bins  -  1 

flag  transferred  from  iread().  1=  end  of  data 

general  purpose  character  used  for  file  names 


file  pointers 


coefficients  for  windowing  waveform  prior  to  FFT. 
read  from  file  "FFTWIN". 

general  counter 

1/df 

temporary  variable  name  of  imaginary  part  of  waveform 
FFT  plus  imaginary  part  of  added  noise  FFT 

imaginary  part  of  the  FFT  of  the  noise,  j  th  noise 
sample,  i  th  frequency  bin.   IT  IS  ASSUMED  THAT  THE 
WARMUP  PERIOD  CONSISTS  ENTIRELY  OF  NOISE. 

counter  used  exclusively  by  iwrite  () 

indicator  transferred  to  actgraf2()  for  graphing. 
l=activity,  0= inactivity 


counter  for  number  of  FFTs  performed 


klug  klug  value  used  to  multiply  variance  to  establish 

threshold 

mag[64]  magnitude  of  the  FFT  of  the  waveform  plus  added  noise 

n  number  of  data   points   in  the   FFT  (128) 

nar  interactively  set  flag  to  indicate  whether  continuous 

display  of  program  variables  k,  varib,  avevar,   thresh 
and  snr  is  desired.     l=desired. 

nisact  counter  for  actual   number  of  declarations  of  activity 

nn  number  of  cycles  indicating  inactivity 

ns  number  of  cycles  indicating  activity 

r[64]  real    part  of  the  waveform  FFT  transferred   from  fft2() 

re  temporary  variable   for  real    part  of  waveform  FFT 

plus  the  real    part  of  the  added  noise  FFT 

rnoise[j][64]  real    part  of  the  FFT  of  the  noise,   j  th  noise  sample, 

i    th   frequency  bin.      IT   IS  ASSUMED  THAT  THE  WARMUP 
PERIOD  CONSISTS   ENTIRELY  OF  NOISE. 

snr  signal-to-noise  ratio,  calculated  for  each  FFT 

ss  sum  of  the  squares  of  the  logs  of  the  magnitudes  of  the 

inner  56  activity  bins,   used  to  calculate  the  variance 

ssn  energy  in   the  outer  4   bins,   used   for  calculating  signal- 

to-noise  ratio 

sss  energy  in   the   inner  56  activity  bins,   used   for  cal- 

culating signal-to-noise  ratio 

sum  sum  of  the  logs  of  the  magnitudes  of  the  inner  56 

activity  bins,   used   to  calculate  the  variance 

suma[4]  the  numerator  in  equation  5),   used  to  calculate  the 

moving,  weighted  mean 

temp  log  base  10  of  the  magnitude  of  the  FFT  bin  currently 

under  consideration. 

thresh  threshold  value   for  declaration  of  activity 

two  interactively  set  flag  to   indicate  whether  2  out  of 

3  criterion   for  activity  declaration   is  desired. 
1 =desired. 

Unext  U         given   in  equation   7)   used   for  calculating  the 

n+1 
moving  variance. 


10 


Uold 

varib 

warmup 

wtave 

x[128] 
y[64] 


U  given  in  equation  7)  used  for  calculating  the 
n 
moving  variance. 

variance  of  the  inner  56  frequency  bins 

interactively  set  number  of  warmup  FFTs 

weighted,  moving  mean  used  to  calculate  moving 
variance 

digitized  waveform  data,  read  from  input  file 

imaginary  part  of  the  waveform  FFT  transferred  from 
fft2() 
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/***********************   PROGRAM :  var diet .  c  ************************ 
* 

*  Written  July,  1984  by  Dr.  J.L  Wayman  to  do  activity/gap  detection 

* 

*  Program  is  written  in  the  C  language. 

* 

*  Subroutines  Required  :  iread(),  f  ft2  ()  ,  actgraf2  ()  ,  iwrite  () 

* 

V 

#include  <stdio.h> 
#include  <math.h> 

#define  PMODE  0644 

short  int  fpr, fpw; 

main  () 


{ 

float 

double 


long  int 


short  int 


int 
char 


r[256], 
Y[256], 
h[128]; 

rnoise[16]  [64] , inoise [16]  [64] , 

mag [64] , suma [4] , 

Uold [4] , Unext [4] , 

sum,  re,  im,  avevar , 

ss , avesnr , 

sss, ssn, snr, 

idf , df , ang, aten, 

denom ,  wt  ave ,  k  1  ug , 

varib, 

temp ,  thresh ; 

two,nisact, 
nn,ns; 

x[128], 

n, 

counter , 

end, warmup, 

i,  j,k, 

countO, count 1, 

count 2, 

nar, 

irec  ; 

isactive; 

fin[20]; 


FILE  *fopen()  ,  *fph; 
FILE  *fopen()  ,  *fp; 


/A************************  OPEN  I/O  FILES  ****************************/ 

if((fp  =  f open ("GRAF", "w"))==NULL) 

{ 

printf ("Error  cannot  open  GRAF\n") ; 

exit  ()  ; 
> 

if((fph  =  fopen("FFTWIN","r"))==NULL) 

{ 

printf ("\nERROR:  cannot  open  %s  for  reading\n", "FFTWIN")  ; 

exit  ()  ; 

y  /*   Read  coefficients  for 

for  (counter =0;  counter<128;  counter ++)  /*  FFT  windowing 

{ 

fscanf  (fph, "%f ", <&h [counter] ) ; 

> 

printf ("\nEnter  output  file  name:"); 
scanf ("%s", fin) ; 

if((fpw  =  creat(fin,PMODE))=-l) 

{ 

printf  ("ERROR,  cannot  open  output  file  for  writing\n")  ; 

exit  ()  ; 
> 

pr  int  f ( " \nEnter  input  f i 1 ename  :  " ) ; 
scanf ("%s", fin) ; 

if  ((fpr=open(fin,0))  =  -1) 

{ 

printf ("\nERROR:  unable  to  open  %s  for  reading\n" , f in) ; 
exit  ()  ;  v 

> 

/********************   INTERACTIVELY  SET  PARAMETERS  ********************/ 

printf  ("\nEnter  added  noise  gain:"); 
scanf  C^f^ang)  ; 

printf ("\nEnter  attenuation  factor  for  moving  averages:"); 
scanf  ('^f'^aten)  ; 

printf ("\nEnter  kluge  factor  for  threshold:"); 
scanf ("%f"/&klug) ; 

printf ("\nEnter  warmup  period  (<=16) : ") ; 
scanf  ("%d" , Awarmup) ; 

printf ("\nNarrative  mode  (l=yes  ,  0=no)  ?") ; 
scanf ("%d"/&nar)  ; 

printf  ("\nTwo  out  of  three  filtering  (l^es,  0=no)  ?")  ; 
scanf ("%d"/&two) ; 


/*************************   INITIALIZE  VARIABLES  *********************/ 

n  =  128; 

irec=0; 

df  =  55.0; 

end  =  O; 

idf  =  1.0/df; 

k  =  O; 

ns=0; 

nn=0; 

nisact=0; 

denom  =l/(l-aten); 

isactive=0; 

countO=0; 

count 1=0 ; 

count 2=0; 

for (i=0;i<4;i++) 

{ 

suma[i]=0; 
Uold[i]=57.0; 
> 

/•A******************  read  AND  PROCESS  WAVEFORM  DATA  *****************/ 

while  (1)  /*  While  there  is  still  raw  */ 

{  /*  waveform  data  to  be  read  */ 

iread (x,  n, &end) ; 

if  (end  =  1)  break;  /*  read  and  window  data     */ 

for  (counter =0; counter <n; counter ++) 

{ 

r  [counter]    =  x  [counter] *h [counter] ; 
y  [counter]    =  O.O; 

fft2(r,y, 128,7,1) ;  /*   Go   to  subroutine   for  FFT   */ 

i  f (k<warmup) 

{ 
for (i=0; i<64; i++) 

{  /*      During  warmup  period  */ 

rnoise[k]  [i]    =  r  [i] ;  /*      store  away  noise    for   later    */ 

inoise  [k]  [i]    =  y[i];  /*      use   in  simulation  */ 

> 
> 

for (i=0; i<64; i++) 

< 

j=k%warmup; 

re  =  r [i]    +   ang*rnoise [j] [i] ;  /*   Add  noise  to   signal    and   */ 

im  =  y[i]    +   ang*inoise[j] [i] ;  /*   take  magnitude  */ 

mag[i]    =  sqrt  ( (double)  (re*re   -t-im*im) )  ; 
> 

sss  =  0;  /*  These  variables  will  be  used   */ 

ssn  =  0;  /*  to  establish  snr  for  this  cycle*/ 


/************************  COMPUTE  ENSEMBLE  NOISE  VARIANCE  *****************/ 

for  (i=0;i<2;i++) 

{ 
temp=loglO(mag[i]) ;  /*   First  two   FFT  bins    */ 

suma [i]    =  aten*suma  [i]    +   temp; 

wtave=  suma [ i ] /denom ; 

Unext[i]    =  Uold[i]    +    (temp* temp   -2 .0*temp*wtave   +  wtave*wtave) ; 

ssn  =  ssn  +  mag[i]  *mag[i]  ; 

temp  =   loglO(mag[i+62]) ;  /*   Last  two   FFT  bins    */ 

suma[i  +  2]    =  aten*suma [i+2]    +   temp; 
wtave  =  suma [i+2] /denom; 

Unext[i+2]    =  Uold[i+2]    +    (temp*temp   -2 .0*temp*wtave   +wtave * wtave) ; 
ssn  =  ssn  +  mag[i+62] *mag[i+62] ; 
> 

/***********  IF  PAST  THE  WARMUP  PERIOD,  COMPUTE  THE  AVERAGE  NOISE  ******** 
************  VARIANCE  TO  THIS  TIME  AND  COMPUTE  THE  ACTIVITY  BAND  ******** 
*************************   VARIANCE   FOR  THI S   CYCLE      ************************ 

avevar  =0; 

i  f  (k>=warmup) 

i 

for (i=0;i<4;i++)  /*  Calculate  average  outer  bin  variance  * 

avevar  =  avevar  +  0.25*(k+l)*  Unext [i] / (denom*k) ; 

thresh  =  (klug  *  avevar)  ;  /*  Set  the  threshold  */ 

isactive=0; 
ss  =  O; 

sum  ?=  0; 

for    (i=4;i<60;i++) 

{ 

sss=  sss   +  mag[i] *mag[i] ;  /*   Calculate  activity  bin  energy  *J 

temp  =   loglO(mag[i] ) ;  /*for   snr   and   find  var   of   log  mag*/ 

ss  =  ss   +   temp* temp; 

sum  =  sum  +   temp; 

> 

varib  =   idf * (ss-sum*sum/ (df+1 .0) ) ; 

if(nar  =  1) 

printf  (M\ncycle%d  varib=%f  avevar=%f  thresh^f^k,  varib,  avevar,  thresh 

/***************  COMPARE  ACTIVITY  BAND  VARIANCE  TO  THE  THRESHOLD  ************ 

if (varib  >  thresh) 

{ 

countO=l; 
ns  =  ns  +  1; 
if  (two=l) 

{ 


if  (countO+countl+count2>=2) 

{ 

isactive=l; 
iwrite (x, 128, &irec)  ; 
> 


> 


else 

{ 

isactive=l ; 

print f ("%d  active\n" , k) 

iwrite (x, 128, &irec) ; 

> 


> 
else 


countO=0; 

nn  =  nn  +  1; 

if(countl  +  count2  =  2) 

{ 

isactive=l .0; 

iwrite (x, 128, &irec) ; 
> 


/*  If  threshold  is  exceeded  and  */ 

/*  2  of  last  3  cycles  also  did  so  */ 

/*  declare  activity  and  write  raw  */ 

/*  waveform  data  to  output  file  */ 


/*  If  threshold  is  not  exceeded  */ 
/*  declare  activity  only  in  the  */ 
/*  event  that  2  of  the  last  3  did*/ 


actgraf2 (fp, isactive) ; 

count2=countl ; 
count l=countO ; 

for  (i=0;  i<4;  i++) 

Uold  [i]  i'  =  aten*Unext  [i]  ; 

k++; 


/*  Dump  activity  status  to  file  */ 
/*  for  printing  later  */ 

/*  update  counters  for  2  out  of  */ 
/*  3  criterion  test  */ 

/*  Attenuate  the  numerator  of   */ 
/*  the  moving,  weighted  variance*/ 

/*  Count  the  number  of  cycles    */ 


/A*******************  CALCULATE  SNR  FOR  THIS  CYCLE  ********************** 
********  SNR  IS  THE  AVERAGE  OVER  EACH  SIX  DECLARATIONS  OF  ACTIVITY  ******/ 

ssn  =  ssn/4; 

sss  =  sss/56-ssn; 

if (isactive==l) 

{ 
if  (nar=l) 

printf ("\n**** 
if  (nisact%6=0) 

avesnr  =  0; 
if  (sss>0) 

{ 
snr  =  10.0*logl0  (sss/ssn) ; 
avesnr  =  avesnr  +  snr/6; 

} 
i  f  (nisact%6=5  &&   nar==l) 

printf ("ave  snr=  %f", avesnr) 


ACTIVITY  ****") 


/*  Calculate  the  average  snr  */ 
/*  over  6  declarations  of  */ 
/*  activity  and  print  out     */ 


nisact++ ; 

> 

/*  Go  back  to  while  statement, */ 

}  /*   go  through  another  cycle    */ 

/* *****************  WHEN  WE'RE  ALL  DONE,  PRINT  THE  NUMBER  **************** 
*******************  OF  CYCLES  AND  COMFIRM  END  INDICATOR  **************** 
*******************  from  FILE  READING  SUBROUTINE  ***************/ 

print f ("k=%d  end=%d\n" , k, end) ; 

> 


DISTRIBUTION  LIST 


DEFENSE  TECHNICAL  INFORMATION  (2) 

CENTER 
CAMERON  STATION 
ALEXANDRIA,  VIRGINIA  22214 


LIBRARY,  Code  0142      (2) 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


RESEARCH  ADMINISTRATION  (1) 
Code  012 

NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


ADJUNCT  PROFESSOR  JAMES  L.  WAYMAN  (15) 

Code  53Ww 

DEPARTMENT  OF  MATHEMATICS 

NAVAL  POSTGRADUATE  SCHOOL 

MONTEREY,  CA  93943 


DEPARTMENT  OF  MATHEMATICS  (1) 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


PROFESSOR  G.  LATTA  (1) 

Code  53Lz 

CHRMN,  DEPARTMENT  OF  MATHEMATICS 

NAVAL  POSTGRADUATE  SCHOOL 

MONTEREY,  CA  93943 


DUDLEY  KNOX  LIBRARY 


3  2768  00336421  7 


